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1.  INTRODUCTION 


A  need  has  long  been  evident  for  a  simple,  fast  numerical  nuclear  fallout 
prediction  model  with  a  capability  approaching  that  of  DELFIC,1*2  hut  with 
greater  accuracy  and  flexibility  than  WSEG-10.3*4  Such  a  code  could  be  used 
for  operational  work,  such  as  damage  assessment  studies,  and  would  be  easily 
usable  by  those  with  less  than  total  commitment  to  fallout  research.  A  pre¬ 
liminary  version  of  a  model,  SIMFIC  (Simplified  fallout  Jnterpretive  Code), 
that  fills  this  need  is  described  here.  The  physical  and  mathematical  bases 
of  the  model  are  described  in  detail  and  validation  results  are  presented. 
Computer  requirements  are  discussed  briefly.  A  complete  description  of  the 
computer  code  is  not  included  since  at  this  time  it  is  in  a  preliminary  stage 
of  development. 
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2.  BASES  OF  THE  SIMFIC  MODEL 


2.1  THE  CENTRAL  PROBLEM 

With  one  notable  exception,  DELFIC,  virtually  all  surface  burst  fallout 
prediction  models  are  patterned  after  the  RAND  stabilized  cloud  model5  which 
was  developed  about  twenty-five  years  ago.  This  type  of  model  is  applicable 
to  explosions  with  energy  yields  greater  than  about  1000  kiloton  (KT)  TNT 
equivalent.  For  such  high  yields  emphasis  is  on  the  cloud  cap  which  contains 
almost  all  of  the  airborne  radioactive  material  at  cloud  stabilization  time. 

On  the  other  hand,  for  low  yield  explosions,  a  few  tens  KT  or  less,  most  of 
the  close-in  fallout  is  in  the  cloud  stem.  When  a  stabilized  cloud  model 
with  an  ill-defined  or  completely  absent  stem  is  used  for  such  cases,  as  has 
all  too  frequently  been  done,  prediction  failure  results  as  is  clearly  shown 
by  recent  studies.6’7 

Experience  with  DELFIC  shows  that  adequate  treatments  of  atmospheric 
transport  and  radioactivity  calculation  are  available.  While  problems  remain 
regarding  particle  size  and  particle  size-activity  distribu  tions,  our  most 
urgent  concern  is  to  define  trajectories  of  particles,  beginning  at  an  early  time 
in  the  fireball,  through  their  rise  along  with  the  cloud,  to  gravity  settle¬ 
ment  and  dispersion  by  the  atmosphere.  This  must  be  done  efficiently  and 
rapidly,  but  with  sufficient  detail  and  fidelity  to  provide  an  accurate  fall¬ 
out  prediction. 

The  key  to  solution  of  this  problem  lies  in  the  well  known  fact  that  the 
altitude  of  a  rising  buoyant  bubble  is  approximately  proportional  to  the  square 
root  of  its  rise  time.  This  has  been  observed,  for  example,  by  Norment  and 
Woolf  for  nuclear  clouds8  and  by  Scorer  for  water  tank  experiments.9  It  was 
used  by  Anderson10  in  an  earlier  fallout  model. 
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In  SIMFIC  the  altitude-square  root  of  time  relationship  is  used  to  set 
up  particle  trajectory  equations  in  the  vertical.  These  are  solved  such  as 
to  develop  equations  from  which  the  maximum  altitude  and  the  time  it  is 
reached  can  be  determined  for  any  particle.  With  this  information  along  with 
ambient  wind  data,  horizontal  trajectory  components  are  added  such  as  to 
define  the  ground  impact  point  for  any  particle  taking  into  account  wind  advec- 
tion  throughout  the  entire  period  of  its  rise  and  settlement.  This  is  done 
for  a  sufficient  number  of  representative  particles  to  adequately  define  the 
fallout  pattern. 

Details  are  presented  in  the  remainder  of  this  chapter. 


2.2  VERTICAL  TRAJECTORY  EQUATIONS 

In  accord  with  the  proportional ity  of  cloud  altitude  with  the  square  root 
of  time,  the  altitudes  of  cloud  (cap)  base  and  top,  z^  and  Zj,  at  time  t  are 
given  by  linear  interpolation  on  /t  between  initial  and  final  altitudes,  zg  , 

ZT,i  and  zB,s  ’  ZT,s’  as 


ZB,i  + 


/t  -  /tT  / 

/r  -  /r  (ZBsS  "  Zb>1,) 


(2.2. IB) 


ZT  =  zT,i  + 


/t  -  /tT 

/r  -  /ft  (Zt>s  ’  Zj'j) 


(2. 2. IT) 


where  the  cloud  rise  is  taken  to  begin  at  time  t..  after  detonation  and  end  at 
time  ts.  Diffe 
top  velocities 


time  t$.  Differentiation  of  equations  (2.2. IB  and  2.2. IT)  give  the  base  and 


(2.2.2B) 


zB,s  '  zB,i 
B  "  2/t(/t;-/t7) 


U,  =  —Lsi - LJ —  (2.2. 2T) 

T  2/t(^  -  /t7) 

Following  the  procedure  used  in  DELFIC,  which  apparently  follows 
Anderson,10  we  assume  a  linear  variation  of  rise  speed  inside  the  cloud 
between  zB  and  Zj.  Thus  the  vertical  velocity  of  an  in-cloud  particle  at 
altitude  z  with  settling  speed  f  is 

df  UB  +  "  ZB^UT  "  UB^ZT  "  ZB^~  f»  z  >  zb»  t  <  ts  (2*2-3) 

Also  following  DELFIC,  we  assume  that  upward  drift  velocity  of  air  below  the 
cloud  decreases  linearly  with  distance  from  the  cloud  base,  which  gives  for 
the  below-cloud  particle  velocity 


dt  =  zlVzB  z  s  zB»  t  <  t$  .  (2.2.4) 

At  this  point  it  is  expedient  to  normalize  the  variables  as  follows: 

?  =  z/(zB,s  "  zB,i)  (2-2-5) 

t  =  /t/(/tJ  -  /t7)  (2.2.6) 

=  dc/dt  =  2ru(/t7  -  ^tT)2/(zB>s  -  zBji), 


v 


(2.2.7) 


and  a  normalized,  average  settling  speed  (see  sec.  2.5)  is  defined  as 


f  =  <  f  >  (/t^-  /t7)2/(zB  s  -  zB  .)  .  (2.2.8) 

In  normalized  form  the  equations  above  become 

=  ^B,i  +  T  ‘  Ti  (2.2. IBn) 

CT  =  ?T,i  +  "  Ti^cT,s  ’  ?T,i^  (2.2.1Tn) 

v  B  =  1  (2.2.2Bn) 

vT  =  cT)S  -  CT>i  .  (2.2.2Tn) 

For  the  in-cloud  particle  we  have 


v  =  1  +  (?  -  ?Bj  -  t  +  t  j )/  (t  -  t.  +  a)  -  2rf;  c  >  cb,t  <  ts  (2.2.3n 

where 


a  '  (?T,i  '  cB,i)/(cT,s  '  ?T,i  '  ’ 

and  for  the  below-cloud  particle  we  have 


(2.2.9) 


v  =  ?/( t  -  t.  +  cB>.)  -  2if ;  ?  <  CB,  t  <  ts 


(2.2.4n) 


After  the  cloud  stops  rising  we  have 


-2rf ;  t  >  t$  . 


(2.2.10) 
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Equations  (2.2.3n),  (2.2.4n)  and  (2.2.10)  can  be  integrated  to  give  the 
vertical  trajectory  equations: 

In-cloud 

C  -  ^  =  (a  +ci  +  CBj)(t  -  T.)/a  -  2f(x  -  t.  +  a)|r  -  r. 

+  (ri  -  a)tn^ - - - ?  >  5g.  r  <  tq,  t  <  (2.2.11) 

Below-cloud 

c  -  E0  “  so(t  ■  To)/(V  Ti  *  5B,i>  ■  2f(t  -  Tf  +  5B,t>  [T  •  To 


+  <ri  '  «B 


C  To  '  T  <  Ts  (2-2 


After  completion  of  cloud  rise 

C  =  ?T=T  -  f(T2  -  x2)  ;  T  >  Ts  .  (2.2.13) 

Here  tq  and  zq  are  the  time  and  altitude  of  separation  of  the  particle  from  the 
cloud  cap.  Next  we  address  the  problem  of  the  determination  of  tq  and  e  . 


2.3  TIME  AND  ALTITUDE  OF  PARTICLE  SEPARATION  FROM  THE  CLOUD 

Fallout  from  the  cloud  cap  occurs  when  c  =  By  substituting  eq. 
(2.2. IBn)  into  the  left  side  of  eq.  (2.2.11)  and  manipulating  the  result 
algebraically  we  arrive  at 

To  ■  Ti  +  a  +  (Ti  '  a)*n(To  “  Ti  +  a)  =  "  tBjl-)/(2af) 


+  (x ,.  -  a)tn(a)  +  a  .  (2.3.1) 


.12) 
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Equations  (2.2.3n),  (2.2. 4n )  and  (2.2.10)  can  be  integrated  to  give  the 
vertical  trajectory  equations: 

In-cloud 

C  -  Cj  =  (a  +ci  +  Cb1)(t  -  xi  )/a  -  2f(x  -  t.  +  a)jr-  -  t. 

+  ( x j  -  a)nn^ - — - -jj|;  ?  >  Sg,  x  <  tq,  t  <  xs  (2.2.11) 

Below-cloud 

C  -  C  =  C  (?  -  X  )/(x  -  X.  +  CR  ,)  -  2f(x  -  T.  +  Cn  -•)  X  -  X 

0  O  0  01  0,1  1  0,11  0 


+  <Ti  -  5b 


?  5  CB’  T0  <  T  C  TS  {2-2-12) 


After  completion  of  cloud  rise 

C  =  CT=T  -  f(x2  -  x2)  ;  x  >  xs  .  (2.2.13) 

Here  xq  and  are  the  time  and  altitude  of  separation  of  the  particle  from  the 
cloud  cap.  Next  we  address  the  problem  of  the  determination  of  xq  and  co- 


2.3  TIME  AND  ALTITUDE  OF  PARTICLE  SEPARATION  FROM  THE  CLOUD 

Fallout  from  the  cloud  cap  occurs  when  c  =  Cg-  By  substituting  eq. 
(2.2. IBn)  into  the  left  side  of  eq.  (2.2.11)  and  manipulating  the  result 
algebraically  we  arrive  at 

tq  -  Ti  +  a  +  (Ti  -  a)tn(xo  -  x.  +  a)  =  (c.  -  CB>i)/(2af) 


+  (x.  -  a)fcn(a)  +  a  .  (2.3.1) 
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This  has  the  form 

£  +  b£n(£)  =  c  (2.3.2) 

which  can  be  solved  for  £  by  Newton  iteration.  Then  tQ  =  £  +  t1-  -  a,  and 
is  found  by  substituting  ro  into  eq.  (2.2. IBn). 


2.4  MAXIMUM  PARTICLE  HEIGHT 

The  quantities  actually  used  in  the  fallout  prediction  are  the  maximum 
altitude  of  the  particle  and  the  time  this  altitude  is  reached. 

For  in-cloud  particles,  set  the  left  side  of  eq.  (2.2.3n)  to  zero,  solve 
for  sm,  the  maximum  altitude,  in  terms  of  x  ,  the  time  at  which  is  reached. 


=  c 


B,i 


-  a  +  2t  f (t 
m  m 


-  Ti  +  a) 


(2.4.1) 


and  substitute  this  into  eq.  (2.2.11).  After  algebraic  manipulation  we  arrive 
at  an  equation  of  the  same  form  as  eq.  (2.3.2),  but  for  which 


£  = 


+  a 


b  =  (t.  -  a )/2 

c  =  (a  +  c.  -  SBji)/(4fa)  +  btn(a)  +  x./2  +  a. 


As  with  eq.  (2.3.2)  this  is  solved  for  £  by  Newton  iteration,  xm  is  recovered 
from  £ ,  and  is  found  by  substitution  of  xm  into  eq.  (2.4.1). 

For  below-cloud  particles,  set  the  left  side  of  eq.  (2.2.4n)  to  zero, 
solve  for  t  in  terms  of  x  . 


2t  f  ( x 
m  m 


-  x. 


(2.4.2) 


substitute  this  into  eq.  (2.2.12),  and  after  algebraic  manipulation  we  again 
arrive  at  an  equation  of  the  form  of  eq.  (2.3.2),  but  with 
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5  =  Tm  ‘  T1  *  5B,i 


b-(.,  -tB>i)/2 

C  =  V[4'<To  '  Tt  *  5B,i>]  *  V2  +  bl"<V  Tf  *  SB,1>  -  Ti  *  cB,i' 


As  before  this  is  solved  for  £  by  Newton  iteration,  etc. 

To  select  appropriate  values  for  x  ,  cm  various  possibilities  must  be 
examined.  (Subscripts  IC  and  BC  denote  in-cloud  and  below-cloud) : 

1.  If  x  >  t  and  t  |Tr  >  x  ,  then  set  t  =  t  and  calculate  c 

osm'lLs  ms  m 

by  substitution  of  xs  into  eq.  (2.2.11). 

2.  If  xo  >  xs  and  xJIC  <  xs,  then  use  (Vcm)ic. 

3.  If  xo  *  xs  and  xs  >  rm|BC  >  ro  and  cJBC  >  cm!rc  or 

Tml  IC  ”  To  use  <V  S.W 

4.  If  xo  <  xs  and  xm|BC  >  tq  and  cJBC  >  cj  IC  or 

Tml IC  >  V  but  Tm!BC  >  Ts  then  set  Tm  =  Ts  and  compute 

by  substitution  of  x$  into  eq.  (2.2.12). 

To  save  computer  time,  zm  and  t  are  computed  only  for  particles  at  the 
base  and  top  of  the  initial  cloud  for  each  particle  size.  For  particles  at 
intermediate  initial  positions,  the  following  interpolation  formula  gives 
acceptable  results 


t  =  t  n  +  k0*85  (t  T  -  t  n ) 
m  m,b  m,T  m^B 


(2.4.3) 


and 
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(2.4.4) 


Computer  plotted  trajectories  for  particles  initially  at  the  initial 
cloud  base  and  top  for  various  particle  sizes  and  weapon  yields  are  presented 
in  Appendix  A. 


2.5  ADVECTION  AND  SETTLING 

A  single  vertical  profile  of  wind  vectors  is  input  to  the  model.  The 
altitudes  at  which  the  vectors  are  defined  are  used  to  stratify  the  atmosphere 
vertically  into  wind  layers,  in  each  of  which  a  unique  wind  vector  applies. 
Advection  and  settling  of  fallout  particles  is  separated  into  two  phases: 

1.  Rise  of  the  particle  to  its  maximum  height  zm  at  time  t  ,  and  2.  settle- 
ment  of  the  particle  from  zm  to  the  ground.  Thus,  if  X  is  the  position  vector, 
relative  to  ground  zero,  of  a  ground  impacted  fallout  parcel  we  have 


where  the  summations  are  over  the  wind  layers  of  depth  Az  between  the  ground 
(z  =  0)  and  the  particle  maximum  height  at  zm.  At  is  the  time  spent  in  a  wind 
layer,  and  ^  is  the  wind  vector  at  height  z. 

During  the  rise  phase,  the  particle  altitude  varies  roughly  as  the  square 
root  of  time.  Therefore  the  time  a  particle  spends  in  the  ith  wind  layer.  At.., 
is  approximately 
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S3* 


At. 

1 


(2.5.2) 


where  a(z2).  =  z?  . -  z ?  and  z.  .  is  the  height  of  the  base  of  the  ith 

i  b,i+l  b ,  i  b  ,i 

wind  layer. 

If  we  define  an  average  settling  speed  between  zm  and  the  ground,  <  f  >, 
and  substitute  eq.  (2.5.2)  into  (2.5.1)  we  obtain 


i  = 


t  z'2 
m  m 


W  a(z" 
z 


)  +  <  f  > 


'E 


W  Az 
z 


z=o 


z=o 


(2.5.3) 


To  avoid  repetitive  calculations  of  the  sums  in  eq.  (2.5.3),  they  are 
precalculated  from  the  ground  to  each  wind  layer  and  stored.  Then  for  each 
particle,  the  partial  summands  that  account  for  transport  to  the  base  of  the 
wind  layer  containing  zm  are  computed,  and  the  balance  of  the  transport  to 
the  ground  is  accounted  for  by  adding  to  them  the  appropriate  presummed 
quantities  that  are  retrieved  from  storage. 

The  average  particle  settling  speed  is  derived  as  follows.  Best11  has 
devised  the  simple  equation 

f  =  f  ebz  (2.5.4) 

0 

for  water  drops,  where  f  is  settling  speed  at  sea  level,  z  is  altitude  (m) 
and 

b  =  2.90  x  10“5;  50  <  6  <  300  pm 
b  =  4.05  x  10“ 5 ;  300  <  6  <  6000  pm 

where  6  is  water  drop  diameter.  This  equation  ulso  gives  adequate  results  for 
fallout  particles. 
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The  average  settling  speed  between  altitudes  7.1  and  z2  is  determined 
by  integration  of  eq.  (2.5.4)  to  give 


f  >  =  fQ(e  2  -  e  *)/  [b(z2  -  zx  )J  . 


(2.5.5) 


In  the  S IMF  I C  code,  a  table  of  f  values  computed  by  the  equations  of 
Beard12  and  Davies13  for  fallout  particles  of  density  2600  kg  m-3  is  stored 
for  75  particle  sizes  (see  sec.  2.8). 


2.6  CLOUD  DEFINITION 

Following  DELFIC,1  the  initial  cloud  is  defined  at  a  time  close  to  the 
fireball  second  temperature  maximum 


t.  =  2.07  W0-13  ,  s  , 


(2.6.1) 


with  radius 


R .  =  108  W°-33  ,  m  , 


(2.6.2) 


with  base  height 


zB,i  =  zGZ  +  zHoB  +  90Wl/3  ’  m  * 


(2.6.3) 


and  with  top  height 


zT  a  =  za  *  +  1.32288  R.  ,  m  . 

1,1  D ,  1  1 


(2.6.4) 


Cloud  stop  and  stabilization  times  are  obtained  by  interpolation  in 
log10  (W)  from  Table  1,  which  is  taken  from  DELFIC  simulations. 


Stabilized  cloud  heights  (m)  are14 


where 


zB,s  *  zGZ  *  zHoB  +  a  “ 


ZT,S  *  zGZ  +  zH0B  tC  “ 


a  =  2228,  b  =  0.3463;  W  <  4.07 

a  =  2661,  b  =  0.2198;  W  >  4.07 

c  =  3597,  d  =  0.2553;  W  <  2.29 

c  =  3170,  d  =  0.4077;  2.29  <  W  <  19 

c  =  6474,  d  *  0. 1650;  W  >  19  . 


TABLE  1 

CLOUD  STOP  AND  STABILIZATION  TIMES 


Yield 

Stop 

Stabil ization 

Time(s) 

Time(s) 

10"3 

300 

421 

lO"2 

300 

421 

10"1 

300 

381 

10° 

300 

382 

101 

300 

422 

102 

280 

663 

103 

200 

783 

104 

160 

787 

105 

150 

991 

(2.6.5) 

(2.6.6) 
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Stabilized  cloud  radius  (m)  is  computed  from 

R$  =  exp(6. 7553  +  0.7381Y  +  0.060308Y2)  (2.6.7) 

where  Y  =  log10(W). 

In  all  cases  above,  W  is  yield  in  kilotons. 

Any  parcel  of  fallout  that  is  inside  the  cloud  cap  at  time  t  is  given  the 

stabilization  radius  Rc.  (The  difference  between  t  and  stabilization  time  is 

s  m 

ignored.)  For  any  size  particle,  the  first  separation  from  the  cap  is  taken  to 
occur  at  time  t.,  and  this  wafer,  which  initially  is  at  the  initial  cloud  base, 
is  given  radius  Ri .  Its  trajectory  is  followed  to  cloud  stabilization,  at 
which  time  its  height,  z  .  ,  which  may  be  below  ground  zero,  is  recorded. 
Parcels  of  fallout  of  this  same  particle  size  at  height  z  are  assigned 
stabilization  radii 


R  ■ R,- +  <z  -  W<RS 


V'<2m 


-  W! 


mm 


<  z  <  z 


B,s 


(2.6.8) 


At  the  initial  time  the  cloud  cap  is  partitioned  into  n  subdivisions  in 
the  vertical  for  each  particle  size  as  illustrated  in  Fig.  1.  n  is  either 
set  to  5,  which  has  been  found  to  be  adequate  for  all  cases  investigated  to 
date,  or 


n  =  10  +  logioW  , 
where  W  is  yield  in  KT. 


(2.6.9) 
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2.7  DISPERSION  BY  AMBIENT  TURBULENCE 

The  same  basic  procedure  is  used  here  as  in  DELFIC.1  Subsequent  to  tm, 
each  parcel  of  fallout  is  taken  to  have  a  Gaussian  distribution  in  the  hori¬ 
zontal  with  variance  given  by15 


o(t)2  =  [°(tm)2/3  +  |  e1/3  (t  -  tm)j  ;  a(tm)  <  o(t)  <  (2.7.1) 

a( t)2  =  1 2(t  -  ^Jfe/a2)1/3  +  3[a(tm)2/aj  ]  ^  "  2|  ’  51  a£ 

(2.7.2) 

where  e  is  turbulence  energy  density  dissipation  rate,  o2  is  the  parcel  vari 
ance  when  its  dispersion  rate  becomes  constant,  taken  to  be  109m2,  and 
a(tm)  =  R/2.  (The  difference  between  tm  and  stabilization  time  is  ignored.) 
Using  Wilkins'approximation  for  e,16  we  define  an  average  e 


and  an  average  settling  time 


-  ‘  ”  ■  V*  f  ' 

where  <  f  >  is  average  settling  speed  (sec.  2.5)  and  zm  is  maximum  parcel 
height  above  ground  (sec.  2.4).  Substituting  these  quantities  into  eqs. 
(2.7.1)  and  (2.7.2),  we  obtain  at  deposition  time,  t^. 


a(td)2  =  [a(tm)2  +  0.26099  zj/3/<  f  >]% 


Z2/3  <  3.83155  <  f  >  JlO3  -  a ( tm ) 2/ 3J 


(2.7.3) 


a(td)2  =  7.8297  x  105  z2/3/<  f  >  +  3  x  106a(tm)2/3  -  2  x  109  ; 


z2/3  >  3.83155  <  f  >  [lO3  -  a(tm)2/3] 


(2.7.4) 


2.8  PARTICLE  SIZE  DISTRIBUTION,  ACTIVITY  CALCULATION  AND  MASS  OF  FALLOUT 

A  table  of  the  75  particle  diameters  larger  than  50  pm  from  a  set  of  100 
computed  by  DELFIC  for  a  lognormal  particle  number  distribution  with  median 
diameter  0.407  and  geometric  standard  deviation  4.0  is  stored.  Corresponding 
to  this  table  is  stored  a  table  of  H  +  1  hour  exposure  rate  activity  fractions 
calculated  by  DELFIC  for  a  0.5  KT  fission  yield  for  fission  types1  U238TN  and 
P239HE  and  averaged. 

A  fallout  prediction  may  be  run  using  the  complete  table  of  75  particle 
sizes,  or  it  can  be  run  using  every  other  entry  (38  particle  sizes),  every 
third  entry  (25  particle  sizes),  or  every  fourth  entry  (19  particle  sizes). 
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K  factors  (Roentgens  m2  hr-1KT-1)  for  seven  fission  types  also  are  stored 
and  any  one  of  these  may  be  selected. 

If  the  initial  cloud  is  divided  vertically  into  n  fallout  parcels  for  each 
particle  size  class  (sec.  2.6),  then  the  area  integrated  activity  at  a  height  of 
three  feet  above  the  ground  from  an  impacted  parcel  of  fallout  in  the  ith  particle 
size  class  is 


Q  =  KWpF.o/n 


(2.8.1) 


where  Wp  is  fission  yield,  F..  is  the  H  +  1  hour  exposure  rate  activity  fraction 
in  the  ith  particle  size  class  and  <t>  is  for  exposure  rate  at  time  t  (hours) 

$  =  t"1-26  , 

and  for  integrated  exposure  (dose)  from  time  tj  to  t2  (hours) 

$  =  (tiO.2*  .  t2°*26)/0.26  . 

Effect  of  height  of  burst  above  ground  zero  is  accounted  for  approximately 
by  multiplying  the  fission  yield  by  the  factor17 


fd  =  (0. 45345)A/6S  ;  A  >  0 


(2.8.2) 


where  A  is  scaled  height  of  burst  in  units  ft  KT-1/3.  This  factor  is  based  on 
a  curve  of  activity  fraction  down  vs.  scaled  height  of  burst  in  Volume  5  of 
DASA  1251. 

The  DELFIC  formulation  of  mass  of  fallout  is  used,1 

m  =  0.07704  W3/3-4(360  +  A) ( 180  -  A)2;  0  s  A  .<  180  (2.8.3) 

where  W  is  total  yield  in  KT  and  A  is  scaled  height  of  burst  in  units 
ft  KT’1/3-4. 

SIMFIC  does  not  have  a  subsurface  burst  capability. 


2.9  DISTRIBUTION  OF  GROUNDED  FALLOUT  PARCELS 

Following  DELFIC1  the  top  and  base  of  each  fallout  parcel  is  transported 
separately  to  ground  impaction.  Then  the  impact  point  coordinates  and  vari¬ 
ances  of  the  two  wafers  are  combined  to  form  parameters  for  a  bivariate 
Gaussian  function  which  is  used  to  distribute  the  parcel  mass  or  activity 
over  the  impact  plane. 

If  the  total  mass  or  area  integrated  activity  in  the  parcel  is  Q,  then  the 
mass  per  unit  area  or  activity  contributed  by  this  parcel  at  point  x,y  in  the 
impact  plane,  u>(x,y),  is 


<o(x,y) 


2tto  a 


(*-xp)2 

2o  2 


(Y~Yp)2~ 

2a  2 


where 


X  =  x  cosa  +  y  sina 


Y  =  y  cosa  -  x  sina 


(2.9.1) 


(2.9.2) 


and  X  and  Y  are  defined  similarly.  The  parcel  central  coordinates,  x  ,  y  , 
P  P  p 

standard  deviations,  o  ,  o  ,  and  orientation  angle  a  are  determined  from  the 

•I 

coordinates  and  standard  deviations  of  the  impacted  parcel  top  and  base  (sub¬ 
scripts  t  and  b  indicate  top  and  base  respectively)  according  to  (see  Fig.  2) 


yp  =  +  V/z 

r  =  /(xt  -  xb)<  +  (yt  -  tb)* 


(2.9.3) 
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Figure  2.  Synthesis  of  a  deposit  increment  standard  deviation  ellipse 
from  the  impact  points  and  standard  deviations  of  the  parcel 
top  and  base  wafers. 


(2.9.4) 


a2  =  (oj.  +  0^  +  r)2/ 4 


a2  =  a, a, 
t  t  b 


cosa  =  I |/r 

sin«  =  (xt  -  xb)(yt  -  yb)/(r|xt  -  xb|) 


(2.9.5) 


To  save  time  in  accumulating  contributions  from  parcels  at  map  points,  a 
table  of  101  exponential  function  values,  spanning  the  range  of  values  ex¬ 
pected  for  the  type  of  map  requested,  is  precomputed  and  stored.  The  expo¬ 
nential  argument  values  of  eq.  (2.9.1)  are  computed  as  indicated,  and  then 
an  approximate  exponential  function  value  is  obtained  directly  from  the  table. 


2.10  FALLOUT  MAP  PREPARATION 

Any  or  all  of  the  fifteen  unique  types  of  fallout  maps  described  in 
Table  2  may  be  computed.  Precisely  the  same  procedure  is  used  by  SIMFIC  as 
is  used  in  DELFIC  for  map  preparation.1 

The  user  may  specify  map  boundaries  and  grid  increments  or  he  may  let 
the  code  supply  these  quantities.  This  is  done  by  the  code  by  considering 
wind  direction,  wind  shear  and  explosion  yield,  and  is  included  for  use  in 
an  exploratory  run  to  obtain  a  first-look  at  a  fallout  map. 
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TABLE  2 


MAP  REQUEST  OPTIONS 


Map  Option  Code  Description 

1  Count  of  fallout  deposit  increments  that  contribute  to 
each  map  ordinate. 

2  Exposure  rate  normalized  to  H  +  1  hour  (Roentgen  hr-1). 

3  Exposure  rate  at  time  H  +  T1  hours,  accounting  for  time 
of  arrival  of  fallout.  (Roentgen  hr-1) 

★ 

4  H  +  1  hour  normalized  exposure  rate  resulting  from 
particles  in  diameter  range  T1  to  T2  micrometers. 
(Roentgen  hr-1) 

5,9  Integrated  exposure  from  H  +  T1  hours  to  infinity, 

accounting  for  time  of  arrival  of  fallout.  (Roentgen) 

6,10  Integrated  exposure  from  H  +  T1  to  H  +  T2  hours,  account¬ 

ing  for  time  of  arrival  of  fallout.  (Roentgen) 

7  Integrated  exposure  from  H  +  T1  to  H  +  T2  hours  assuming 
all  fallout  has  arrived  by  H  +  T1  hours.  (Roentgen) 

8  Integrated  exposure  from  H  +  T1  hours  to  infinity  assum¬ 
ing  all  fallout  has  arrived  by  H  +  T1  hours.  (Roentgen) 

11  Mass  of  fallout  per  unit  area  (kg  m-3). 

12  Mass  of  fallout  per  unit  area  deposited  from  H  +  T1  to 

H  +  T2  hours  (kg  m-3 ). 

13  Mass  of  fallout  per  unit  area  deposited  by  particles  in 
diameter  range  T1  to  T2  micrometers,  (kg  m-3 ) 

15  Time  of  onset  of  fallout,  (s) 

16  Time  of  cession  of  fallout,  (s) 

17  Diameter  of  smallest  particle  deposited,  (pm) 

18  Diameter  of  largest  particle  deposited,  (pm) 


n 

A  "normalized"  calculation  is  one  in  which  it  is  assumed  that  all  fallout  is 
deposited  by  H  +  t  regardless  of  actual  deposition  time. 


3.  COMPUTATION  REQUIREMENTS 


SIMFIC  is  a  FORTRAN  code  presently  operational  on  the  CDC  6600  computer. 
Including  extensive  glossaries  and  comments,  the  code  uses  1545  statement 
cards,  which  is  about  3/4  of  a  box  of  cards.  It  uses  less  than  16,000  central 
memory  storage  words  and  no  peripheral  storage. 

To  obtain  a  benchmark  computing  time,  the  five  predictions  discussed  in 
the  next  section  were  done  in  one  run  using  five  vertical  cloud  subdivisions 
(sec.  2.6)  and  nineteen  particle  sizes  (sec.  2.8).  Computing  time  was  10.5 
seconds.  The  same  five  predictions  were  done  in  one  run  using  the  SEER  code.1 
Twenty-five  particle  sizes  were  used  and  the  code  set  the  number  of  vertical 
cloud  subdivisions  (n  =  3  +  sinW);  identical  maps  were  produced  in  terms  of 
boundaries,  grid  increments  and  numbers  of  points.  Computing  time  was  7.4 
seconds.  Thus,  SIMFIC  is  almost  as  fast  as  SEER,  but  produces  much  superior 
predictions.6’7  (It  also  uses  about  one-third  as  much  central  processor 
storage  as  SEER.)  For  damage  assessment  work,  where  many  predictions  are 
required  for  fixed  explosion  yield,  substantial  savings  of  computing  time 
could  be  realized  without  sacrifice  of  prediction  accuracy  and  with  only 
modest  increase  in  central  processor  storage  usage. 

Data  requirements  of  the  code  are  minimal:  yields  (total  and  fission), 
height  of  burst  above  ground  zero,  ground  zero  height  above  sea  level,  fission 
type  (optional),  data  for  a  single  wind  profile,  map  type  code  and  quantities 
T1  and  T2  if  required  (Table  2),  map  boundaries  and  grid  intervals  (optional*) 
and  some  simple  control  information. 


★ 

The  code  has  a  capability  to  set  map  boundaries  and  grid  intervals  based  on 
wind  direction,  wind  shear  and  yield. 


4.  VALIDATION 


4.1  DISCUSSION 

Predictions  by  SIMFIC  and  DELFIC  are  compared  with  observed  H  +  1  hour 
normalized*  exposure  rate  maps  for  the  five  test  shots  described  in  Table  3. 
DELFIC  predictions  were  executed  as  discussed  in  reference  7,  using  the  data 
listed  there.  SIMFIC  predictions  were  similarly  executed  using  H  hour  winds 
only,  with  five  cloud  subdivisions  (sec.  2.6)  and  nineteen  particle  sizes 
(sec.  2.8).  Observed  fallout  patterns  were  taken  from  DASA  1251. 

Three  methods  of  comparison  of  fallout  patterns  are  used: 

1.  Visual  comparison  of  contour  maps. 

2.  Comparison  of  contour  areas,  and  hotline  lengths  and  azimuths.** 

3.  The  Rowland-Thompson  Figure-of-Merit  (FM).19  (Appendix  B) 

These  are  roughly  in  order  of  importance. 

Statistical  data  are  in  Table  4  and  the  contour  plots  are  on  pp.  30 
through  44.  The  contours  were  drawn  by  a  30-inch  Calcomp  plotter,  and 
members  of  each  observed-predicted  set  are  to  the  same  scale. 

TABLE  3 

TEST  SHOT  DATA 


Shot 

Total 

Yield 

(KT) 

Fission 

Yield 

(KT) 

HOB 

(m) 

Altitude 
of  GZ 
(m) 

Site 

Johnie  Boy 

0.5 

0.5 

-0.584 

1570.6 

NTS+ 

Jangle-S 

1.2 

1.2 

1.067 

1284.7 

NTS 

Small  Boy 

low 

- 

3.048 

938.2 

NTS 

Koon 

150. 

- 

4.145 

0.0 

Bikini 

Zuni 

3380. 

- 

2.743 

0.0 

Bikini 

+Nevada  Test  Site 

A  "normalized"  exposure  rate  map  is  constructed  on  the  assumption  that  all  local 
fallout  is  down  at  the  specified  time,  regardless  of  its  actual  deposition  time. 

*Hotline  length  is  defined  as  the  furthest  distance  from  ground  zero  on  a  contour, 
and  hotline  azimuth  is  the  angle,  measured  clockwise  from  north,  to  the  point  of 
furthest  distance  from  ground  zero  on  a  contour. 


TABLE  4 


COMPARISON  OF  OBSERVED  AND  PREDICTED  FALLOUT  PATTERN  STATISTICS 
Observed/DELFIC/SIMFIC 


FM 

DELFIC 

Contour 

Hotl  ine 

Test  Shot 

SIMFIC 

(Roentqen  hr'1) 

Area (km2 ) 

Lenqth(km) 

Johnie  Boy 

0.182 

1000 

0.278/0.029/0.018 

1.38/0.32/0.32 

0.152 

100 

0.539/0.774/0.781 

2.73/2.58/2.69 

50 

1.271/1.787/1.828 

4.10/4.13/5.14 

58(42)/61  (44  )* 

28(3)734(13)* 

Jangle-S 

0.483 

500 

0.117/0.144/0.156 

0.69/1.00/1.60 

0.437 

300 

0.386/0.316/0.435 

1.50/1.23/2.23 

100 

1.437/2.242/1.278 

3.74/5.87/3.40 

35 

3.114/5.077/4.093 

5.06/7.68/9.37 

40(45)/22(18) 

43(42)769(48) 

Sma 1 1  Boy 

0.308 

1000 

0.216/0.047/0.144 

1.00/0.25/0.60 

0.533 

500 

0.528/0.135/0.375 

1.62/0.56/1.15 

200 

0.942/0.564/1.031 

2.22/1.69/2.17 

100 

3.75/1.10/2.26 

5.66/3.72/4.02 

50 

9.03/4.38/7.59 

8.10/6.47/8.17 

63(59)725(24) 

44(36)/20(15) 

Koon 

0.287 

500 

32.0/26.0/44.0 

10.2/12.5/14.9 

0.325 

250 

122/87.3/116 

17.3/24.2/24.1 

100 

550/261/374 

41.0/39.5/41.6 

33(40)725(18) 

22(22)/29(20) 

Zuni 

0.105 

150 

474/2239/2854 

98/78/77 

0.189 

100 

2761/3619/4566 

125/96/97 

50 

6187/6660/9365 

138/121/142 

30 

10950/9913/15190 

177/153/180 

105(16)7164(52) 

17(16)/12(9) 

f 

Mean  absolute  percent  errors:  DELFIC/SIMFIC.  The  values  in  parentheses  are  calculated 
including  the  data  for  the  highest  activity  level  contours.  See  footnote  next  page. 


Aziwuth(deg) 

359/  0/  0 
345/344/344 
343/343/344 


342/353/351 
346/354/351 
1/355/352 
6/355/  6 


71/  66/  80 
73/  80/  75 
72/  73/  71 
72/  74/  68 
75/  72/  66 


18/  0/  0 
15/  4/  0 

17/  3/  1 


12/337/354 

17/337/356 

27/338/355 

33/340/359 


without 


Prediction  accuracy  is  seen  to  be  good,  particularly  for  the  low  yield 
shots.  Overall  mean  absolute  percent  errors*  for  contour  area  and  hotline 
length  are: 

Contour  Area  Hotline  Length 
DELFIC  61(42)  32(26) 

SIMFIC  59(31)  32(21) 


The  quantities  in  parentheses  are  computed  with  the  data  for  the  highest 
activity  level  contours  excluded.  The  highest  level  contours  are  particu¬ 
larly  difficult  to  predict,  usually  being  in  the  region  affected  by  throwout 
and  induced  activity  in  and  around  the  crater.  SIMFIC  and  DELFIC  do  not 
address  this  portion  of  the  activity  field  since  fallout  is  a  negligible 
contributor  to  casualties  there.  SIMFIC  prediction  accuracy  is  insensitive 
to  numbers  above  five  of  vertical  cloud  subdivisions.  Increasing  the  number 
of  particle  size  classes  smooths  the  patterns  and  eliminates  the  multiple 
contour  closures  in  the  Jangle-S  prediction,  but  appears  not  to  warrant  the 
increased  cost  in  computing  time. 

For  these  cases  SIMFIC  is  shown  to  be  at  least  equal  to  DELFIC  in  pre¬ 
diction  accuracy.  It  is  important  to  emphasize  that  this  level  of  competency 
in  SIMFIC  has  been  achieved  without  a  posteriori  adjustment  or  calibration  of 
any  aspect  of  the  model  so  as  to  improve  agreement  with  any  observed  fallout 
pattern. 

The  three  low  yield  shots  were  executed  at  the  Nevada  Test  Site,  and 
their  fallout  patterns  were  measured  over  land.  For  this  reason,  observed 
patterns  for  these  shots,  though  not  highly  accurate,  may  be  considered  to 

★ 

For  n  observed-predicted  data  pairs,  mean  absolute  percent  error  is 


; 


be  superior  to  the  patterns  of  the  high  yield  shots  which  were  executed  on 
Bikini  Atoll  in  the  South  Pacific.  Not  only  are  the  fallout  fields  of  the 
high  yield  shots  very  large,  which  adds  to  measurement  problems,  but  most 
of  the  fallout  from  these  shots  fell  into  water.  Even  so,  most  of  the  Koon 
pattern  area  was  covered  by  an  array  of  fallout  collection  stations,  so  this 
pattern  is  probably  reasonably  accurate.  Zuni,  on  the  other  hand,  is  a 
special  case.  The  falTout  pattern  used  here  is  exclusively  downwind  of  the 
atoll  and  was  determined  by  an  oceanographic  survey  method  that  was  known 
to  be  inaccurate.  The  close-in  pattern  in  the  region  of  the  atoll  is 
available,  but  contains  no  closed  contours  so  could  not  be  used  here;  thus 
the  high-activity  portion  of  the  observed  pattern  for  this  shot  is  ignored, 
and  this  alone  must  account  for  a  substantial  portion  of  the  disagreement 
between  observation  and  prediction  for  this  shot,  particularly  with  regard 
to  contour  areas  and  contour  overlap  (Table  4).  In  addition,  we  have  the 
following  problem. 

Predictions  for  these  high  yield  shots  are  expected  to  be  inferior  to 
those  for  these  low  yield  shots.  This  is  because  both  of  the  high  yield 
shots  were  detonated  over  coral  soil,  and  in  the  case  of  Zuni,  a  large  but 
uncertain  amount  of  sea  water  was  lifted  by  the  cloud.  The  particle  size 
distribution  used  for  these  predictions  is  typical  of  fallout  produced 
from  the  siliceous  soil  found  at  the  Nevada  Test  Site.  We  have  not  suc¬ 
ceeded  in  developing  a  distribution  appropriate  for  coral  and  coral-sea 
water  mixtures. 

More  details  concerning  the  prediction  calculations  and  test  shot  char¬ 
acteristics  are  in  reference  7. 


29 


OBSERVED 
SMALL.  BOY 


REFERENCES 


1.  H.  G.  Norment,  "DELFIC:  Department  of  Defense  Fallout  Prediction  System. 
Volume  I  -  Fundamentals,"  Atmospheric  Science  Associates,  DNA  5159-1 

(26  October  1979). 

2.  H.  G.  Norment,  "DELFIC:  Department  of  Defense  Fallout  Prediction  System. 
Volume  II  -  User's  Manual,"  Atmospheric  Science  Associates,  DNA  5159-2 
(26  October  1979). 

3.  G.  E.  Pugh  and  R.  J.  Galiano,  "An  Analytical  Model  for  Close-In  Opera¬ 
tional-Type  Studies,"  Weapons  Systems  Evaluation  Group,  WSEG  RM  No.  10 
(15  October  1959).  AD  261  752. 

4.  R.  B.  Mason,  "Description  of  Mathematics  for  the  Single  Integrated  Damage 
Analysis  Capability  (SIDAC),"  National  Military  Command  System  Support 
Center,  NMCSSC  TM  15-73  (1  July  1973).  AD  913  164L. 

5.  W.  W.  Kellogg,  R.  R.  Rapp  and  S.  M.  Greenfield,  "Close-In  Fallout," 

J.  Meteor.  14,  1  (1957). 

6.  H.  G.  Norment,  "Analysis  and  Comparison  of  Fallout  Prediction  Models," 
Atmospheric  Science  Associates,  unpublished. 

7.  H.  G.  Norment,  "Evaluation  of  Three  Fallout  Prediction  Models:  DELFIC, 
SEER  and  WSEG-10,"  Atmospheric  Science  Associates,  DNA  5285F 

(16  June  1978). 

8.  H.  G.  Norment  and  W.  Woolf,  "Studies  of  Nuclear  Cloud  Rise  and  Growth," 
Technical  Operations,  Inc.,  unpublished. 

9.  R.  S.  Scorer,  "Experiments  on  Convection  of  Isolated  Masses  of  Buoyant 
F1 uid,"  J.  Fluid  Mech.  2,  583  (1957). 

10.  A.  D.  Anderson,  "A  Theory  for  Close-In  Fallout  from  Land-Surface  Nuclear 
Bursts,"  J.  Meteor.  18,  431  (1961). 

11.  A.  C.  Best,  "Empirical  Formulae  for  the  Terminal  Velocity  of  Water  Props 
Falling  Through  the  Atmosphere,"  Quart.  J.  Roy.  Meteor.  Soc.  76,  302 
(1950. 

12.  K.  V.  Beard,  "Terminal  Velocity  and  Shape  of  Cloud  and  Precipitation 
Drops  Aloft,"  J.  Atm.  Sci.  33,  851  (1976). 

-  * —  V  -  ■  - - 

i-'HLCLuUG  r>.,J £  bWi (-WOT  IljJ&D  W 


47 


REFERENCES,  continued 


13.  C.  N.  Davies,  "Definitive  Equations  for  the  Fluid  Resistance  of  Spheres," 
Proc.  Phys.  Soc.  (London)  Jtf,  259  (1945). 

14.  E.  F.  Wilsey  and  C.  Crisco,  "An  Improved  Method  to  Predict  Nuclear  Cloud 
Heights,"  Ballistics  Research  Laboratory,  unpublished. 

15.  J.  J.  Walton,  "Scale  Dependent  Diffusion,"  J.  Appl.  Meteor.  12,  547 
(1973). 

16.  E.  M.  Wilkins,  "Decay  Rates  for  Turbulent  Energy  Throughout  the  Atmos¬ 
phere,"  J.  Atm.  Sci.  20,  473  (1963). 

17.  R.  L.  Showers,  "Improvements  to  the  PROFET  Fallout  Prediction  Program," 

Ballistics  Research  Laboratories,  BRL-MR-2095  (February  1971).  AD-883  280. 

18.  H.  Lee,  P.  W.  Wong  and  S.  L.  Brown,  "SEER  II:  A  New  Damage  Assessment 
Fallout  Model,"  Stanford  Research  Institute,  DNA  3008F  (May  1972). 

AD- 754  144. 

19.  R.  H.  Rowland  and  J.  H.  Thompson,  "A  Method  for  Comparing  Fallout  Pat¬ 
terns,"  DAS I AC,  G.E.  -  Tempo,  DNA  2919F  (April  1972). 


APPENDIX  A 


PARTICLE  TRAJECTORIES 


Computer  drawn  trajectories  in  the  vertical  computed  by  the  methods  of 
secs.  2.1  -  2.4  are  presented  for  various  particle  sizes  and  weapon  yields. 
All  particles  initially  are  at  the  initial  cloud  base  or  top.  Heights  and 
times  are  normalized  according  to  eqs.  (2.2.5)  and  (2.2.6).  Yields  are  in 
KT  and  particle  diameters  in  pm.  Horizontal  and  vertical  lines  are  drawn 
across  the  graphs  to  intersect  at  points  of  separation  from  the  cloud  cap 
and  maximum  height. 
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APPENDIX  B 


FALLOUT  PATTERN  COMPARISON  BY  THE  FIGURE-OF-MERIT  METHOD 


Rowland  and  Thompson19  developed  this  method  for  comparison  of  pairs 
of  fallout  contour  maps  by  computation  of  a  single  index,  the  FM,  that  is 
a  measure  of  contour  overlap  between  them.  For  each  contour  common  to  the 
patterns,  the  area  overlapped  and  the  area  not  overlapped  is  calculated. 

The  areas  are  weighted  by  the  average  radiation  level  between  successive 
contours.  Sums  over  all  contours  of  weighted  overlapped  areas  and  weighted 
total  areas  are  computed,  and  the  FM  is  the  ratio  of  the  two  sums.  For  com¬ 
pletely  overlapped,  perfectly  matched  patterns,  FM  =  1;  for  no  overlap,  FM  =  0. 

Mathematically,  FM  is 


FM 


y(ri +  ri-i> 
w  2 


N 


<ri  +  ri-i) 


i»l 


0>i  -  «M> 

(S(  -  A,.!) 


where 

N  is  the  number  of  contours  in  the  patterns.  The  summations 
are  from  highest  contour  to  lowest 


r..  is  activity  of  the  i  th  contour  ( Roentgen/ hr ) ,  rQ  =  10  rj 

a ^  is  common  (i.e.,  overlapped)  area  for  the  i  th  contours. 

a  =  0. 
o 

Ai  is  total  area  of  the  i  th  contour.  The  summation  in  the  de¬ 
nominator  is  computed  for  both  patterns,  and  the  largest  sum 
is  used.  A  =  0. 


The  FM  has  been  found  to  have  limited  utility  as  a  measure  of  fallout 
prediction  accuracy.  This  is  mainly  for  two  reasons.  First,  and  most 
important,  is  that  being  a  measure  of  overlap,  the  FM  is  strongly  biased 
in  favor  of  overprediction;  that  is,  it  favors  predictions  that  cover  a 
large  area,  and  therefore  overlap  the  observed  pattern,  regardless  of 
other  considerations.  Second,  the  FM  method  imposes  no  penality  for  miss¬ 
ing  or  extra  contours;  contours  not  common  to  both  patterns  are  simply 
ignored. 
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